ROBUST ERROR ESTIMATES FOR STABILIZED FINITE ELEMENT 

APPROXIMATIONS OF THE TWO DIMENSIONAL 
NAVIER-STOKES EQUATIONS WITH APPLICATION TO IMPLICIT 
LARGE EDDY SIMULATION 

ERIK BURMAN* 

Abstract. Wc consider error estimates in weak parametrised norms for stabilized finite element 
approximations of the two-dimensional Navier-Stokes' equations. These weak norms can be related 
to the norms of certain filtered quantities, where the parameter of the norm, relates to the filter 
width. Under the assumption of the existence of a certain decomposition of the solution, into large 
eddies and fine scale fluctuations, the constants of the estimates are proven to be independent of both 
the Reynolds number and the Sobolev norm of the exact solution. Instead they exhibit exponential 
growth with a coefficient proportional to the maximum gradient of the large eddies. The error 
estimates are on a posteriori form, but using Sobolev injections valid on finite element spaces and 
the properties of the stabilization operators the residuals may be upper bounded uniformly, leading 
to robust a priori error estimates. 
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1. Introduction. In this paper we will be interested in stabilized finite element 
methods in the context of so called implicit large eddy simulation (ILES), see [2]. This 
is a numerical approach to the computation of turbulent flow where no modelling of 
the Reynolds stresses is performed on the continuous level. Instead the Navier-Stokes' 
equations are approximated numerically using a method that dissipates sufficient en- 
ergy on the scale of the mesh size. This eliminates the buildup of energy that creates 
spurious oscillations in any energy conservative approximation method. It has been 
argued that the truncation error of such methods by itself may act as a subgrid model 
[TJ [18] and there exists numerical evidence that ILES methods work for the simulation 
of two dimensional turbulence, provided back scatter effects are not strong [T6J. There 
is also numerical evidence of the potential for adaptive LES/DNS driven by adaptive, 
stabilized finite clement simulations, see [12l[T3j and [21) . 

Our objective in this paper is to provide a numerical analysis for stabilized finite 
element methods under minimal regularity assumptions and to provide sufficient con- 
ditions on the exact solution for the derivation of rigorous error estimates that are 
independent of both the Reynolds number and Sobolev norms of the exact solution. It 
is well known that provided the exact solution is sufficiently smooth the approximate 
solution Uh of the Navier-Stokes' equations on velocity-pressure form can be proved 
to satisfy estimates of the type 

||«-«fc||L»(n) <e llVulU hi\u\ L 2 (I . H 2 m , (1.1) 

if a consistent stabilized finite element method with piecewise affine approximation is 
used. Here u denotes the flow velocity and h the mesh size. See [lit fT] for examples of 
analyses of Navier-Stokes' equations on velocity-pressure form and [15l|T7] for analyses 
on vclocity-vorticity form. Wc also give a proof of (jl.ll) for one the methods proposed 
herein in appendix. Here we use the notation a < b for a < Cb with C a constant 
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independent of the physical parameters of the problem, unless they can be expected 
to have 0(1) contribution, it can also include some dependence on initial data, that 
may be assumed to be O(l). We will also use a ~ b for a < b and b < a. Note that 
there is no explicit dependence on the viscosity in the estimate For this estimate 

to be useful the included Sobolcv norms must be small, which rarely is the case in 
the high Reynolds number regime and hence the dependence of the viscosity enters 
in an implicit manner. The purpose of the present paper is to propose an alternative 
approach, where the estimate is indeed independent of the Reynolds number, both in 
the sense that the estimate is free from inverse power of the viscosity in the upper 
bound, but also that the dependence on unknown Sobolev norms of the exact solution 
is strongly reduced. It does not seem possible to eliminate this dependence completely, 
due to the possible presence of backscattcr. Observe that the presence of ||Vu||oo in 
the exponential of (jl.ip reflects the effect of diverging characteristics in the transport 
equation and is present already in the linear convection-diffusion equation at high 
Pcclct numbers. We may define a timescale for the flow separation, r := ||Vu||~ , . 
The reason this time scale becomes so small is that it will be the smallest timescale 
of the flow and hence equal to the micro time, because of fine scale fluctuations of 
the velocity. From the physical point of view it is argued that LES will be successful 
for flows where both the quantities of interest and the rate- controlling processes are 
determined by the resolved large scales, see Pope [20] for a discussion. We will use this 
idea as a starting point for our assumptions on the flow. 

To derive error estimates for a numerical method we need the following: 

- continuous dependence on data, independent of the exact solution; 

- some smooth quantity that we can apply approximation estimates to. 

At a first glance both these prerequisites appear to fail for the two-dimensional Navier- 
Stokes' equation. The first fails because of the presence of the exponential factor and 
the second fails because Sobolev norms of the exact solution can be huge for small 
viscosities. The following three points allow us to break this deadlock: 

1. the use of a parametrized weak norm, corresponding to measuring the error 
in filtered quantities of the solution; 

2. introduction of an assumption on the structure of the exact solution that is 
sufficient for an implicit large eddy simulation to be robust; 

3. a stabilized finite element method, giving enhanced a priori control of residual 
quantities in the high Reynolds regime. 

The idea of measuring error in filtered quantities was considered in [9j [10] , but the 
estimates were not robust in the Reynolds number and the constant included high 
order Sobolev norms of the exact solution. In [5] weak norm estimates were used 
in order to derive robust estimates for the Burgers' equation, where the constant in 
the right hand side only depended on initial data. The second point, which was not 
necessary in the case of the Burgers' equation, reflects the difficulty to characterise 
the solution structure in higher dimension and the ensuing need of some structural 
assumptions in order to rule out strong backscattcr effects. The third point allows us 
to control residual quantities independent of the viscosity. 

As a first approximation it is reasonable to assume that for a solution to be 
amenable to large eddy simulation, there are relatively smooth eddies, with large 
associated Reynolds number, containing the bulk of the energy and small scale fluctu- 
ations that may vary rapidly in space, but carries a negligible part of the energy. To 
make this precise, we assume that there exists a decomposition of the exact solution 
in the spatially slowly varying large scales and an arbitrarily rough fine scale, with 
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small energy. 

u = u + u', u,u' G W 1,00 (Q). 

We then assume that the Reynolds number associated to the large scales, Re may be 
large, but ||u||-j/yi,°<>(Q) ~ 1, whereas for the fine scales ||w'||via,«>(q) may be large, but 
the energy small. To give a precise meaning to small here, we introduce a global time 
scale for the flow, defined using the large scales 

t f ■= IMIwA.c^q) ~ i. 

This is in agreement with the statement that rate controlling processes are deter- 
mined by the large scale. Using the viscosity coefficient we may make the following 
assumption on the energy content of the small scales 

h'Wlc ~ v/t F - 
The length scale based on \u'\ and Tp writes 

l':= KlTF ~ Kl™ 



and it follows that the small scale Reynolds number is 

Re > := *L „ „ L (L2) 
v \u'\v 

Alternatively one may assume that the fine scale Reynolds number is one and that 
the large scale characteristic time, is the globally relevant time scale and then derive 
the bound on the energy. We will refer to the above as the large eddy assumption. 
Under this assumption we prove the following bound on the approximate velocities 

sup \\u~u h \\ L 2 {n) <h^. (1.3) 
te(o,T) 

The hidden constant in the above estimate only depends on initial data (maximum 
initial vorticity) and the mesh geometry, but is independent of the Reynolds number 
and Sobolev norms of the exact solution. The discussion is limited to two space 
dimensions and hence we do not properly speaking address the question of turbulent 
flows. 

Let us end this introductory discussion by emphasising that what we compute is an 
approximation to the solution of the Navier-Stokes' equations. For this approximate 
solution we can prove that provided Tp is not too small, corresponding to slowly 
varying large scale velocity field, the filtered part of the vorticity is stable under 
perturbations resulting in robust error estimates in weak norms for vorticity. Using 
these estimates we may then control the i 2 -norm of the velocity error as shown above. 
Herein our main concern will be the high mesh Reynolds number case 

Re h := > 1, 

v 

where Uq := 0)||i=c(n) ~ 1 denotes a characteristic velocity of the flow, but 

many results are independent of the mesh Reynolds number. It will always be ex- 
plicitly stated when a result only holds in the high Reynolds regime. If the local 
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Reynolds number is low, other approaches than those presented herein might be more 
appropriate. Let us also point out that another feature of our estimates is that they 
provide the first error estimates with an order in h for nonlinear stabilization schemes, 
satisfying a discrete maximum principle, in two space dimensions. 

We will consider the Navier-Stokes' equations written on vorticity-velocity form. 
Let £1 be the unit square and assume that the boundary conditions are periodic in 
both cartesian directions. The L 2 -scalar product over some space-time domain will 
be denoted (•, -)x with associated norm || • \\x where the subscript may be dropped for 
X = f2. Define the time interval / := (0,T) and the space-time domain Q := SI x /. 
The equations then writes, 

dtOJ + V-(itw) — vAuj = 0, in Q, 

-A* = uj in Q, (1.4) 
u = rot "3/ in Q, 
w(x,0) = uj Q , 

with ujq G L°°(il). The associated weak formulation takes the form, for t > 0, find 
(w,*) S ff^Q) x ij!(0) ni t (0) such that 

(d t uj,v)n + a(u;cj,v) = 0, (1.5) 

(V*,V$)n = (w J *)n, (1-6) 
u = rot 

for all (v, $) G iJ 1 (Sl) x if 1 (J7) n i*(Sl), where the semi-linear form a(-; •, •) is defined 

by 

a(u; u), v) := (V-(vu), v)n + (vVu, Vu)n- 

This problem is known to be well-posed, but a priori error estimates on the solution 
are in general strongly dependent on the viscosity coefficient reflecting the possible 
poor stability of the equations in the high Reynolds number regime. 

2. Finite element discretization. Let {Th]h>o be a family of affine, simplicial 
Delaunay meshes of f2. We assume that the meshes are kept fixed in time and that the 
family {Th}h>o is quasi-uniform. Mesh faces are collected in the set T . For a smooth 
enough function v that is possibly double- valued at F G J- with F = dT~ n dT + , 
we define its jump at F as [u] :~ v|t- ~ w It+? an d we fix the unit normal vector 
to F, denoted by up, as pointing from T~ to T + . The arbitrariness in the sign of 
{vj is irrelevant in what follows. Define Vh to be the standard space of piecewise 
polynomial, continuous periodic functions, 

Vfi := {v h G H l (Q) : v h \ K G P k (K);VK G T h \ v h periodic in x and y}. 

The set of gradients of functions in Vt will be denoted by 

Wl' 1 := {w h = Vv h ; v h G V*}. 

Let L := L 2 (fl) and set £» := {q G L; J n q = 0}. Let V* ■= V£ Ci L*. We let tt l 
denote the L 2 -projection on V£ and 7ry the H 1 -projection 



(Vttvu, Vv h ) n = (Vm, Vv h ) n Vv h G Vh and / (n v u - u) dx = 



n 
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We recall that the following approximation estimates hold for ir^ and 7Ty, 

\\tt l u - u h \\ + h\\V(TT L u -u)\\ < c h s \u\ s , with 1 < s < k + 1 (2.1) 

and 

\\tt v u - u h \\ + h\\V(TT V u -u)\\ < dh s \u\ s , with 1 < s < k + 1. (2.2) 

We consider continuous finite elements with k = 1 to discretize the vorticity uj in space 
and k = 1, 2 for the stream function The discrete velocity is given elementwise by 
uh\k '■= rot vf/j := (dy^h, —d x ^h)- Note that using this definition V • Uh = in fi, 
i.e. the discrete velocity is globally divergence free. We discretise in space using a 
stabilized finite element method. For t > find {uh^h) £ V h l x Vj , with I = 1,2, 
such that 

(d t oj h , v h ) M + a(u h ; u h , v h ) + s(u h ; oj h , v h ) = 0, (2.3) 
(V* fc , V$,0o - (w h> $ ?i )o = 0, (2.4) 
Uh - rot = 0, 

for all (vh 1 §h) G Vft x K and with initial data wo ■= 7Tl^>(-, 0). Here s(-; •, •) denotes a 
stabilization operator that is linear in its last argument and (dtUJh,Vh)M denotes the 
bilinear form defining the mass matrix, this operator either coincides with (-, -)q or is 
defined as (■, -)n approximated using nodal quadrature, i.e. so called mass lumping. 
We will assume the stabilization term satisfies 

inf \\h^{u h ■ \7io h - v h )\\ < s(u h ,uj h ]uj h )^ < hh (U + \\u h \\L<*>(n))\\VcJh\\ (2.5) 

vhGVh 

and 

s{u h ,uj h ;v h ) < h?(U a + \\uh\\L^(U))s(uh,^h;^h)^\\^v h \\. (2.6) 

The formulation (|2.3p - (|2.4p satisfies the following stability estimates 
Lemma 2.1. 

sup + 2||i/*Vwi,||^ + 2 f s(u h ;u h ,w h ) dt < \\u H {; 0)\\ 2 M , (2.7) 

tei J i 

and if exact integration is used for (•, -)u> 

\\u h (-,T)\\lf + 2\\viu, h \\l = \\u h (-MM-2 [ s(u h ;uj h ,* h ) dt, (2.8) 

Jo 

K(-,*)||L-(n) <cJwh(-,t)||L, ( n), 9>2 (2.9) 

and for 1 = 1, 

/ \\Vd t u h \\dt< / (h-^(U + \\u h \\ L oo m )s{u h ,u} h ,u} h )^ +vh- 2 \\Vuj h \\)dt. (2.10) 
Jo Jo 

Proof. Inequality (|2.7I) is immediate by taking Vh = in (|2.3[) . Inequality (|2.8[) 
is obtained by taking vt = &h m the equation (|2.3I) and deriving the equation p. 41) 
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in time and taking = uj^. For the inequality (|2.9[) . consider the auxiliary problem, 
— = LUh in Vt and note that by [35] there holds 

||«fc(-,«)|U-(n) £ ll*(-^)ll^.«(n) 

and adapting the analysis of |19j we have for the (simpler) case or periodic boundary 
conditions, 

\\^{-,t)\\ W i,^ {u) < c q \\u h (-,t)\\ Lq{n) , q>2. 

To prove (|2.10[) finally we introduce a function £ Vh such that 

(€h,Vh)M = (^d t u h ,Vv h )a, Vv h £ V h , 

it follows by taking Vh = £,h and using the Cauchy-Schwarz inequality followed by an 
inverse inequality that 

(Zh^hVM-MhWZh^WdtVcOhW. (2.11) 

Observe that by norm equivalence on discrete spaces the L 2 -norm defined using nodal 
quadrature is equivalent to the consistent L -norm. Taking Vh = £h hi (|2.3[) yields 

\\dtVuj h \\ 2 = -(u h ■ Vuj h ,^ h )n - (vVoJh, V^)n - s(u h ,u> h ,£ h ). 

We may then apply the Cauchy-Schwarz inequality in the second term of the right 
hand side and (|2.6j) in the last term, followed by inverse inequalities on ||V^|| and 
the estimate (|2.11[) . For the first term we write, using the properties of ^ and the 
bound 

IK,a)M-K,a)n| < (h 2 \Vv h \,\V£ h \) Q , 

\{u h ■ VoJh,€h)n, I < \{uh ■ Vwf, - v h ,^ h ) n \ + \(d t Vu h ,Vv h )n\ + (h 2 \Vv h \, \V£ h \)n 

Since both Uh and Vwj are constant per element Vi^ \k = V(u/, — Uh ■ Vw/, ) | k ■ Using 
inverse inequalities and the bound (|2.11[) on ^ we have 

\{u h -Vw h -v h ,£ h )\ + \(d t Vu h ,Vv h )\ + (h 2 \Vv h \,\V£ h \) < h^WdtVLJhWWuh-Vwh-Vhl 

The claim follows by the inequality (|2.5|) . (|2.1ip and finally by integrating in time. □ 
It follows from (|2.8[) that the method is energy consistent if s(uh', ujh, ^h) = 0. Taking 
the difference of the formulations (|1 .5[) - (| 1 . 6j) (with v = vt) and (|2.3[) - (|2.4[) and 
setting e u = u) — u)h and = — "f^, the following consistency relation holds 

(d t e u +u-Ve w + rot eyVuj h ,v h ) n + {vWe u , Vv h ) u = (d t uj h ,v h ) M - {d t uj h ,v h )n 

+ s(u hl uj h ;v h ) in Q, (2.12) 
(Ve*, V$fc) n - (eo,, = in Q. 

As mentioned in the introduction, if the solution (u,u>) is smooth one may prove an 
error estimate that is robust with respect to v using standard linear theory and per- 
turbation arguments. For the methods we consider herein, this result is an extension 
of the works in |17j and [7] and we state it here only with the dominant terms present. 
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For the readers convenience we briefly outline the proof using one stabilization oper- 
ator (defined in equation (|5.5[) ) in the appendix. 

PROPOSITION 2.2. Let (u,u) be a smooth solution of (|f .5j) - (|1.6j) and (u h ,uj h ) 
be the solution of (|2.3[) - (|2.4[) . where the stabilization operator satisfies the additional 
weak consistency property 

s(u h ;n L uj,Tr L uj) < c(u,uj)h^ 

then for I = 1 , 2 

||(u - u h )(-,T)\\ + || (w - w h )(-,T)\\ < c u (hi \u)\ L * {I , mm + h l \^\ L ^ {I . H ,+i {n)) ) 

where c w := e" Vlj; " LOC(< " T . In addition there holds for the stabilization operator 

s(u h ;u] h ,uj h )^ < c w (/i^H L 2 (/ . ff 2 (0)) + h l \ty\ LX{I . H i+i {n)) ). 

Observe that the exponential factor here depends on ||Vc«;||ioo(o)i compared to 
|| Vu||l°°(o) in (jl.ll) . This is the prize we pay for estimating the L 2 -error of the 
vorticity. As we shall see below, the use of weaker norms for the estimation of u)h 
allows us to revert back to the exponential factor of and under the large eddy 
assumption, the exponential growth is moderate. 

3. Dual problem. From the consistency relation (|2.12[) we deduce the following 
(homogeneous) perturbation formulation for the evolution of (e u ,e$) 

(d t e w + u-Ve w + rot e^-Vw/j, ipi)a + (iA7e w , Vpi)n = in Q, (3.1) 

(Ve*, V(p 2 )n - {e u ,y>2)a. = in Q, (3.2) 

where ip±, if2 are the solutions to a dual adjoint perturbation equation related to the 
continuous equation (| 1 . 5[) - (| 1 and the discretization (|2.3p ~ (|2.4p . Since the jump of 
the tangential derivative of u)h is zero, we may integrate by parts in (|3.1|) . to arrive 



at the dual adjoint problem 

—dtf\ — u-Vipi — ip2 — vAipi = in Q, (3-3) 

— A(/?2 — Vcj/j-rot <pi = in Q, (3.4) 

<pi(x,T) = £o(x) in (3.5) 



where £,o(x) is some initial data to be fixed later, the choice of £o determines the 
quantity of interest. 

A key result for the present analysis is the following stability estimate for the dual 
adjoint solution 

Proposition 3.1. The following stability estimate holds for the solution (ipi, ip^) 

of (33) - i|33]t . 

sup ||V^i(-,t)|| + ||^£>V||q < e^||V£o|| (3.6) 
tei 

J \\VM-M dt<e% J \\u} h \\ L ao (a) #||V£ || (3.7) 
where Tp is defined in the proof. If the large eddy assumption holds Tp ~ 1 . 
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Proof. First multiply (|3.3[) by — A(pi and (|3.4[) by tpi and integrate over Q* := 
fix (t*,T), where t* is a time to be chosen. By summing the two relations we obtain 

(d t ipi, Atp^Q. + ju-Vipi , Atpi) Q , + (Vuh-rot i pi,tpi) Q ^ + ||z/5Av?i||q, = 0. 
h h h ~~u~ 

We will now treat the terms term by term. First note that by integration by 

parts first in space and then integration in time we have 

h = -\[, Jt H v ^M)H 2d * = \\\V<Pi{;t*)\\ 2 ~ |||Ve || a . 

The second term is handled using the decomposition of u in the large scale and fine 
scale component and then an integration by parts only in the large scale part. Here 
Vsit denotes the symmetric part of the gradient of the vector u. 

h = -((Vsu- -(V-u)l2x2)Vpi, V<Pi)q* - (u'-V^i, A(^i)q« 



< 



(A(?2,u',^)V^i) T -V^i dxdt + -\\viA<pi\\ 2 Q ,, 



where A(u, u', v) is a two by two, symmetric matrix defined by, 

1 _ It 
A(«, u',v) = -V s u+ -V-uI 2X 2 + -z-u u. 

We now define the global timcscale tf of the flow by 

Here er+ denotes the largest positive eigenvalue of the matrix. This results in a 
nontrivial minimization problem in L°° . We leave the precise study of this problem 
for further work and here simply observe that by computing the eigenvalues of the 
symmetric part of the gradient tensor we may write 

(tf)' 1 < MJ(u,u') 

u 

where 

J(u,u') := sup (^IKd^Si - d X2 u 2 ) 2 + (d x . 2 ui + <9 Xi u 2 ) 2 ||l°°(q) + IKIIi°°(n)) • 

We observe that the global stability does not depend on the divergence component 
or the rotational of u, only on the other two components of the velocity gradient 
matrix. Since u' = u — u, it follows that we can minimize over all large scale vector 
fields u <G [VF 1,00 (fi)] and the infimum value obtained is the optimal timescale of the 
flow. Under the assumptions made in the introduction, that ||u||wi-o°(n) ~ 1 an d 
v _1 Il u 'lli^(f2) ~ 1; f° r a ^ we immediately deduce that tf ~ 1. 

By an integration by parts and by using the relations V-rot ip = and V</?-rot ip = 
we have 



I 3 = -(w/jV-rot <^i,^i)q» - {uj h rot <p 1} V^i)q» = 0. 
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Collecting the results for I± — I 3 we have 

HV^iM*)|| 2 L + W^A^Wl, < r^HV^II 2 ^ + HVfcHi. 

The inequality for ipi follows after a Gronwall's inequality and by taking the supremum 
over t* , resulting in 



sup ||V^(-,t)|| 2 + M^VIIq < e*||V£|| s 



Elliptic regularity has been used for the second term. 

For the bound on ip2 multiply equation (|3.4j) by if2 and integrate over f2, 



\\V<p 2 {-,t)\\ 2 = -(w h rot <pi,V(pa)a < |M-,t)||L~ (n) ||V^i(-,t)||||V<p 2 (-,t) 
Then divide by ||Vy>2 (•,£)!!, integrate in time and use that 



Wh(-,i)IU»(n)||V¥>i(-,i)|| dt< / \\io h (;t)\\ L ^ m dt sup ||V^(-,t)||. 
i Ji tei 

Finally use equation (|3.6j) to bound the term in ||Vy>i(-, t)\\. □ 

Note the dependence of ujh in the bound (|3.7|) . This appearance of a finite element 
function in the stability estimate shows that the global stability depends on the mono- 
tonicity of the approximation scheme. However as we shall see, strict monotonicity is 
not necessary, only L°°-control of the vorticity. 

4. A posteriori and a priori error estimates for the abstract method. 

Let e u = lu — Uh and let the filtered error e w be defined as the solution to the problem 

- 5 2 A£ U + i u = e u . (4.1) 

We introduce a norm on e u such that |||e w |||^ := H^Ve^H 2 + |e w || 2 = (e^e^^. This 
norm coincides with the L 2 -norm for <5 = and is related to the 7J _1 -norm for 5 = 1. 
By choosing 5 := 5(h), i.e. by reducing the filter width with the mesh size, we obtain 
a family of norms that become stronger as the mesh size is reduced. 

Using the above norm and the relations (|2.12[) , p.3[) - (|3.4j) as well as the stability 
result of Proposition ^ . 1 1 we may derive a posteriori estimates for the filtered quantity 
iuj- We here derive the result for the abstract finite element element method p.3p - ()2.4l) 
and then show how these estimates can be transformed into a priori error estimates, 
depending on the properties of the stabilization operator s(uh,^h',Vh). The use of 
weak norms and stabilized finite clement methods in the following estimates draws on 
ideas from [14] and OH]. 

Theorem 4.1. (A posteriori error estimates) 



i=0 



with 

TZq := || (w - Uh)(- 



TZi := / inf \\h^ (u h -VuJh - v h ) || dt, 

Jo v h £V h 
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TZ 2 := nm(/i > v5T*)||i/*[n F .Vwfc]||.Fxi, 



^3:= / ||wh(-,t)||i<»(n) dtmin(c sup |j*/ l (-,t)|| A ,o,ci/i 2 sup 
Jo tei tei 

where 

||*fc(-,t)lk« := \\h s ln F -\7y h (-,t)]y+ mf yi [ ]T ||fc*+'(Att h (., t) - v h )f K 



Vh 



\KeT h 



3 f T 

Ki := ft* I ||<9 t Vwh|| dt 
Jo 



and 

TZ 5 := (U + ||u/i||ioo(Q)) / s(u h ;u> h ,u> h )2 dt. 

Jo 

The term IZ4 is omitted if the consistent mass matrix is used. For the velocities we 
have the estimate, for all t G I , 

\\(u-u h ){;t)\\ < (\\* h (-,t)\\ A!h + \\\(Q-Q h )(-,m\i) (4-3) 

where the second term in the right hand side may be a posteriori bounded by taking 
5=1 in (|4~2|) . 

Proof. By the definition of e w we have, taking £0 = e w in (|3.5[) . 



\\\eu\\\l = (e w (T), <pi(T))n + (e w , -d t tpi - u-Vipi - vA<pi)q 
+ (e*, -A(p 2 - Voj h ■ rot ipi) Q 
= (e w (0), <pi(0))n + (d t e w + u-Ve w + rot e*-Vo; fe , + (VVe„, V<^i)n 
+ (Ve*, Vtp 2 )n - (e w , ^ 2 )n- 
Using now the consistency relation (|2.12p we obtain 

|||e w ]]|a = (e w (0),(^i -ttl<Pi)(',0))(1 + (d t e + u-Ve + rot e*-Vw/,,^i -itl<Pi)q 
+ (>Ve, V(^i - 7r L ^i))Q + (Ve*, V(<^2 - n<^ 2 )) Q - (e, y 2 - n^ 2 ) Q 
- {d t uJh,K L Lpi) M .Q + (d t Uh,ir L <pi)Q - s(uh,Uh\nL<Pi), 



where II : M> will be taken as cither ttl or 7ry. Using the equations (jl.5 

(|1.6[) and the dchnitions of the projections 7T£ and 7ry we deduce for II := Try, 

We^Hl! = (e w (0), (<^i - 7Ti,^i)(-,0))n - (u h -Vu> h - v h ,<pi - tvlVi)q 
- (i/Vw/,, V(y>i - ttl<Pi)q + {uh,<P2 - ^V^P2)Q 

-{dtUJ h: -K L Lp- l ) Mt Q + {dtU) h ,'KL^l)Q- / s{uh,w h ;-K L ipx) dt, 



and similarly for II := 7Tl, 

lllewllli = (e w (-,0), (<pi - n L <pi)(;0))n - {u h -Vu h - v h ),Vi - ttlvOq 

- (vVWft, V((y5l - TTlV>i)q - (Wfc, V(<^ 2 - KL<P2))q 
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After some standard manipulation including integrations by parts, Cauchy-Schwarz 
inequalities, trace inequalities the approximation results (|2.1[) and (|2.2I) we may con- 
clude, for II := 7ry, 



W\eJ\ 2 s <(^j (\K(;0)\\ + j Q v MJhHu h -Vuj h -v h )\\dt 

+ mm{h^T^)\\^[n F -yLo h j\\jr XI + ah^ sup \\u h (; t)\\ / \\u h (-,t)\\ L ^ m 

tei Jo 

+ \\d t Vu h \\ dt+ (U + \\uh\\L°°(Q)) / s(u h ;uj h ,uj h )^ dtj 



dt 



x ( S up||5V^i(-,t)|| + ||^ J DVilk)- 
tei 

If II := ttl the fourth term on the right hand side is replaced using 



(V*fc,V(¥>2 -7t z ,^ 2 ))q < ( ) cosup||* fc (i)||A,o / \\5Vtp 2 (;t)\\ dt, 



tei 



followed by the bound ([3. 7(1 on f2- The estimate (|4.2[) now follows by taking the 
minimum of the two expressions and noting that by (|3.6|) 



sup ||5V¥?iC-,t)|| + ||<5i^L>V||q ^ e-^ |||e w ||| 5 . 

The velocity estimate (|4.3[) is obtained by noting that, with := \& — 

\\u - UhW 2 := ||Ve^|| 2 = (Ve$, V(e* - 7r L e*)) + (e w ,7Tie*). 
Using the equation (|1.6[) we have 

|Ve*|[ 2 = (V^h, V(e^ - 7TLe*)) + (w, e$) - (wh, 7rz,e*) 

= (V* ft , V(e* - 7r L e*)) + (w - w h , e*) 

Let e be the solution of (|4.ip with 5=1. Then 

||u - "hll 2 = (V*/,, V(e* - Tr L e^))n + (Ve w , Ve*)o + (e w ,e*)n- 

By an integration by parts in the first term, followed by a Cauchy-Schwarz inequality 
and the Poincarc-Fricdrichs inequality in the last term we may write 

ll« - u h \\ 2 < |l4lV*h]||^||/»~»(e* - n L e*)\\r 

+ ( \\h(AV h - v h )f K ) \\h-\e 9 - 7r ie *)|| 

\K£T h ) 

+ 111(a) - - u h )||. 

By elementwise trace inequalities and the approximation property (|2.ip we have 
||/i~i(e* - 7r L e*)||jr + ||ft _1 (e* - 7r L e*)|| < \\u - u h \\ 
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by which we conclude. □ 

If the stability properties of the stabilized method arc sufficient, these a posteriori 
error estimates translate into a priori error estimates. We propose two strategies for 
this. One using stability concepts based on Sobolev injections for discrete spaces and 
one based on monotonicity, applicable to monotone stabilized finite clement methods 
and monotone implicit large eddy methods. The advantage of the former is that it 
allows the derivation of a priori error estimates for quasi linear terms s(u^; w/ l5 v/i) and 
the use of the consistent mass matrix. The latter technique on the other hand allows 
for the derivation of a priori error estimates with precise control of the constants in the 
estimates. We will use the notion of the discrete maximum principle (DMP) and the 
associated, DMP-property of the forms defining a finite element method introduced 
in 0. 

Proposition 4.2. Assume that the mass (•, -)m is evaluated exactly and that in 
addition to (|2.5[) and (|2.6[) the following stability estimate holds for all t > 0, 

||w/,||i=o ( n) <c(h)(\\u) h \\ + s(u h ,Lj h ,u) h )i). (4.4) 
Then there holds for all e > 0, 

111(0 - &h)(T)\\\s < (jjs) 2 (ceh-* + c(h)h^) 

and 

\\(u-u h ){;T)\\ < inf \\(u-v h )(;T)\\+e^h?(c e h-c + c(h)h*). 



Proof First we recall that |K(-,0)|| < ||w(-,0)||. Then by (J231 and (|2~7| 

/ inf \\ht(u h -Vu h -v h )\\ dt<T?( / a{u h ;w h ,u h ) dt) 2 < T*||w(-,0)||. 
Jo v h£V h \J ) 

Using an elementwise trace inequality and (|2.7|) we also have 



min(h,v?Ti)\\viln F -Vu) h }\\F X i < hv\\v$Vw H \\ Q < fc*||w(-,0)||. 
For 72-3 we use the discrete Sobolev injection ([4.4)1 to deduce 



h* sup \\u h (-,t) 
tel 



\\uh(-,t)\\ L ™(n) dt 



<hnw h (;0)\\<h) / (|K|| 

Jo 



+ s(u h ,LU h ,LL> h )2) dt 

<hh(h)\\uj h (;0)\\ 2 . 



The only remaining term is the stabilization term, which is not innocent since we do 
not have an a priori bound on the factor ||nfe||i,°°(Q). Here wc use (|2.9I) to deduce, for 
all i > and q > 2, 

||«fc|U°°(fi) J s(u h ;u) h ,u) h )2 dt <c q \\u h \\ Lq{n) T^yj s(u h ;uj h ,uj h ) dtj 
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and by a global inverse inequality and the bound (|2.7|) we may conclude 




and the estimate follows taking e = (q — 2)/q. Note that the constant c q explodes as 

The bound in the L 2 -norm for the velocities follows as before from the vorticity 
estimate using, with Ch denoting the Clement interpolant, 

||u - u h \\ 2 = ||Ve*|| 2 = (Vej, V(tf - C h *)) - (e u , V h - C h V) 

< \\u - u/,||||V(* - C h V)\\ - (-Ae u + e u , (C h V - * h )) 

< ||u - «fc||||V(* - C h *)\\ + \\\u - - *h||Hi(n)- 

We conclude by using the H 1 -stability of the Clement interpolant, a Poincarc inequal- 
ity and finally by dividing both sides with ||u — Uh\\- □ 

Proposition 4.3. (A priori error estimate using monotonicity) Assume that 
Ren > 1, that the mass (-,-)m is evaluated using nodal quadrature, that the form 
a(uh]uih,Vh) + s{uh]U)h,Vh) has the DMP property as defined in JSj/ and that (|2.5jl - 
(|2.6p are satisfied as well as the assumptions of Lemma \2.1\ Then there holds 

|||(C-<5 fc )(T)|||,<e^ ^ 

and 

\\(u-u h )(;T)\\ < inf \\(u-v h )(-,T)\\ + e%h$. 

Proof. The terms TZq — 72-2 are bounded as in the proof of Proposition 14.21 Since 
by assumption the spatial discretization of (|2.3|) has the DMP property and the mass- 
matrix is evaluated using nodal quadrature, we know from [U [6] that 

IM|l~(q) = ||wfc(-,0)||i«. ( n). 

Hence by (|2.9|) ||«/i|[l>»(q) < Coo|[wh(-, 0)||ioo(nj. We may then use these L°°-bounds 
together with the stabilities of Lemma 12.11 to upper bound the remaining residual 
quantities of (|4.2I) . Using (|2.9[) and (|2.7[) we immediately have 

fc*sup||w fc (-,t)|| f ||w/i(-,t)||icc (n) d*</iir||w h (.,0)||||w fc (.,0)||iao ( n) 
tei Jo 

For the residual term resulting from the mass lumping we have using the stability 

(USD 

hi J \\d t Vu h \\ dt < T5(C/ + \\u h \\ L oo (Q) )(J {s{u h -u h ,Lo h ) + ||^ Vw/l || 2 ) dt) 2 

<t--u!m-MV 
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The remaining contribution from the stabilization is bounded as before using the 
maximum principle and (|2.7|) . The proof of the L 2 -norm estimate on the velocities is 
identical to that of Proposition 14.21 □ 

Note that only the proof of Proposition 14.31 uses the assumption Reh > 1 and only to 
control the non-consistent mass term. This constraint is likely to vanish if the method 
is analysed using techniques appropriate for parabolic problems, since mass lumping 
is known to be stable for dominant diffusion (see for instance [24]). 

5. Stabilized finite element methods. The estimates of Theorem 14.11 holds 
for any finite element method on the form (|2.3|) - (|2.4p . Indeed by taking s(-; ■, •) = 
the standard Galerkin method is included. This means that in general the effect of 
stabilization can be observed only a posteriori, by observing smaller residuals for the 
stabilized formulations. In Propositions 14.31 and 14.21 we propose a priori estimates 
derived from the a posteriori error estimates under special assumptions on the prop- 
erties of the stabilizing terms. These can be proven to hold only for stabilized finite 
element methods, since the standard Galerkin method does not allow for a control 
of the second term of the right hand side of (|4.2[) independently of the viscosity, nor 
can (|2.5p and (|4.4p be proven to hold. In this section we will suggest some stabiliza- 
tion operators that satisfy the assumptions necessary for the results of the abstract 
analysis to hold. We will consider the following cases: 

1. linear artificial viscosity, in which the numerical viscosity is increased so that 
the mesh Reynolds number always is one. Using a lumped mass matrix 
together with anisotropic viscosity we may design the scheme to satisfy a 
discrete maximum principle, giving a priori control of ||w/i||.l«>(q)- When the 
consistent mass matrix is used one may proved that (|4.4[) holds giving once 
again a priori estimates, at the price of a logarithmic factor. 

2. high order stabilization, we propose to stabilize the jump of the streamline 
derivative. This scheme does not yield a maximum principle, so the residuals 
can not be completely a priori bounded. The scheme has some interesting 
conservation properties for two-dimensional Navier-Stokes' computations that 
we will point out. If a nonlinear stabilization term is added and mass-lumping 
is used the solution may be made monotone and the a priori error estimate 
of Proposition ^. 31 holds, this time with the possibility of higher order conver- 
gence in the smooth portion of the flow. Finally if the consistent mass matrix 
is used and stabilization is added also in the crosswind direction, an estimate 
of the type (|4.4j) can be shown to hold leading to a priori error bounds using 
Proposition (|4.2jl . 

5.1. Methods using consistent mass matrix. We consider first the stabi- 
lization method obtained by penalizing the jumps of the streamline derivative over 
element faces. We use the exact mass matrix in (|2.3j) and the stabilizing operator 

s sd (u hl uj hl v h ) := 7 U Q 1 (h 2 F lu h -Vw h j, [uh-Vv h j) F - (5.1) 

FG.F 

For this formulation the following stability estimates hold 
Lemma 5.1. 

sup||^(.,i)|| 2 +2||^V^||2 +2 7 cV 1 ||M U/l -V W/! .]||^< IM-,0)|| 2 (5.2) 
tei 

and if the consistent mass matrix is used, 

jM.,T)|| 2 + 2||^|| 2 Q HK(.,0)|| 2 (5-3) 
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Proof, the proof of (|5.2I) is an immediate consequence of (|2.7|) and the definition 
(|5.1[) . The inequality (|5.3[) follows by observing that 

s ad (u h ,u h ,^ h ) = 7 XI (^^fK-VwJ, [rot * h -Vtf/J|)ir = 0. 
fg.f 

□ 

Observe that the method dissipates enstrophy but conserves energy exactly as the 
physics of the problem suggests. Using known results on interpolation between dis- 
crete spaces it is also straightforward to show (see [8] ) , 

inf \\h?(uh • Vwd - v h )\\ 2 < s sd (uh,Uh,Uh)- 

Vh&V h 

Unfortunately this stabilization operator can not be shown to satisfy (|4.4p . For this 
we need the stabilization to act also in the crosswind direction. We therefore propose 
the following two stabilization operators, the first is the standard artificial viscosity 
method 

s av (u h ;tL> h ,LL> h ) := (jhCUo + \u h \) 2 U ' 1 \7u) h ,Vv h ) (5.4) 

and the second is a modification of (|5.1[) where also the crosswind gradient is penalized 
defined by 

Scd{uh,Uh,Vh) ■= Ssd(uh,Uh,Vh) +7i y2 U Q h^ K / \n F ■ VwJ^f • Vv/J ds (5.5) 

Ken JdK 

Observe that the first part of s cd ensures the satisfaction of (|2.5|) and as we shall see 
the second part is necessary for (|4.4[) to hold. 

PROPOSITION 5.2. Both stabilization operators (|5.4|) and (|5.5|l satisfy (|2.5|l and 
(|2.6|) . T/ie stabilization operator s av (-; •, •) satisfy (|4.4|) with c(/i) ~ /i _ 2(l -f |?o<?(/i)|) 

and Scd(-; ■, •) satisfy (|4.4j) wit/i c(/i) ~ h 4 (1 + |Zog(/i)|), // > 0. 

Proof. The proofs of (|2.5j) - (|2.6|) are consequences of the Cauchy-Schwarz in- 
equality and in the case of s c d trace inequalities. To prove (|4.4[) we note that in two 
space dimensions there holds (see [25] ). 

||wh||i«.( n) < (1+ |log(/i)|)||o;h|| H i (a) . 

This allows us to conclude for s„„. For s r ri we use that 



Vw h ||<(^ / |[Vw h -nji-llKI ds 



A Cauchy-Schwarz inequality followed by a trace inequality in the right hand side 
leads to 

i 

HVwfcH < h~^\\oj h \\K\\h^{nf-Wuh\\\aK\ 2 < h-^(\\u h \\+s cd (u h ;u; h ,u} h )i). 

Ken 

□ 

Since the assumptions of Proposition 14.21 arc satisfied, we may conclude that the 
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method (|2.3|) - (|2.4I) using the stabilization (|5.4j) statisfy the a priori error bounds for 
e > 

UK* - "h)(T)\\\ S < (J^j 2 (c e h-< + 1 + | Iog(ft)|) 

and 

i 

ll(« - u h )(-,T)\\ < inf - v h )(-,T)\\ + ( |r ) ' (c e /i" e + 1 + | \og(h)\). 

Similarly we have the following estimates if the stabilization ()5.5|) is used. 

i 

O h )(T)|||i < (*L\ 2 ( Ce h-< + (1 + | log^Q^) 

and 

||(« - u h )(-,T)\\ < inf ||(u - v h )(-,T)\\ + e^h?(c e h-* + (1 + | log(/i)|)/r ^). 

We see that if we take \i = 1 in (|5.5[) we get the same order for the two methods, 
however if we want the method to have optimal convergence for smooth solutions we 
choose fi = 2 and I = 2, resulting in an a priori convergence order of 0(/ji) in the 
non-smooth case. 

5.2. Monotone methods. Since the consistent mass matrix is non-monotone 
we herein only consider methods using lumped mass. Monotone methods can also be 
designed using a nonlinear switch that changes the local quadrature as a function of 
the solution u>h so that the consistent mass is used away from local extrema to reduce 
the dispersion error known to haunt mass-lumping schemes, such methods are beyond 
the scope of the present paper. 

5.2.1. Linear artificial viscosity. A monotone method using linear artificial 
viscosity is obtained by taking (see [3]) 

s{u h ,uJh,Vh) ■= r y^{aiax(U , \\u h \\ L c* iK) )h 2 K } (Vu/, x n F , Vv h x n F ) F . (5.6) 

K FedK 



Then the estimates (|2.7|) - (|2.8[) hold and we observe that there exists positive constants 
ci , C2 such that 

cilllufcl^Vw/jH^ < s{u h ,u h ,u h ) < C2\\\uh\^h^WuJh\\h- 

Let the mass matrix be evaluated using nodal quadrature so that the matrix cor- 
responding to (•, -)a/ is diagonal. We may use the theory of [51 [B] to prove that 
the operator a(wh, U/J + s(uh, ujf, , v^) has the DMP-property and hence the following 
discrete maximum principle holds 

IKIIl~(q) = ||w h (-,0)|| £a =(n). 

This requires the parameter 7 to be chosen large enough, however it does not require 
any additional acute condition on the mesh, since the discretization of the Laplace op- 
erator results in an M-matrix on Delaunay meshes. Since by the maximum principle, 
H«fclU°°(Q) < ll w / l (- I 0)|| L oo (n) we have 

IIKI^Va^lH < KUl^ / s(u h ,oj h ,u) h ) dt < \\u> h (;0)\\M (5-7) 
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which proves (|2.5p with Vh — 0. It is straightforward to prove also (|2.6[) . Comparing 
with Proposition 14.31 we conclude that the assumptions are satisfied and hence that 
the Proposition holds for (|2.3p - (|2.4[) with stabilization given by (|5.6[) and the mass 
matrix evaluated using nodal quadrature. 

5.2.2. Nonlinear artificial viscosity. Here we assume that I = 1 so that both 
uih and are discretized using piecewise affine elements. We propose a stabilization 
term consisting of one linear part and one nonlinear part. The role of the nonlinear 
part is to ensure that the form a(-; •, •) + s(-; •, •) has the DMP property. The linear 
part is necessary to ensure that the inequality (|2.5|) holds. We define 

s(u h ;u) h ,Vh) ■= s sd {u h ;uj h ,v h ) (5.8) 
+ 72^^ 2 ^ R F (u hl uj h )(sigii(\7uj h x n F ) 1 Vv h x n F ) F (5.9) 

K F £dK 

where 

R F {u h ,uj h ) := ||ith||i«.( Aji .)(l + ^c7 1 |I u ?i||l~(a f ))™-f([«. j f • Vw/,]) 
with A F := \J KeTh . KnF ^K and 

m F ({n F ■ VwftJ) = max \\[n F ■ VcjJUf'- 
F'er 
F 'edK';K'n F=F 

It is shown in [5] that with this definition a(ui l ',coh,'Vh) + s(uh',^h,Vh) has the DMP- 
property for 72 large enough. Since the bounds (|2.5[) - (|2.10p also hold, the assumptions 
of Proposition 14.31 are satisfied and its estimates hold. We conclude that for the 
methods defined by mass lumping and the stabilization operators (|5.6p or (|5.8p the 
following estimates hold 

1 

\\\{6j-Q h )(T)\\\s<e% 

and 

||(u-u h )(.,T)||< inf ||(«- Wfc )(.,r)||+ e ^ 

6. Conclusion. We have shown that under a certain structural assumption on 
the solution of the two dimensional Navier-Stokes' equation one may derive robust 
error estimates with an order in h, independent of both the Reynolds number and 
high order Sobolcv norms of the exact solution. Robustness is obtained for a class 
of stabilized finite element methods. The estimates are both on a posteriori form, 
and on a priori form, providing an upper bound on the error. Due to the strong 
assumptions on the mesh the present a posteriori error estimates are not immediately 
suitable for use in adaptive algorithms, but a more detailed analysis may allow the 
mesh assumptions to be relaxed. If the solution is smooth we also prove that optimal 
convergence may be obtained, provided the stabilization operator is weakly consistent 
to the right order. 

Observe that it is natural that the LES estimate has much poorer convergence 
order, since we may assume no smoothness of the exact solution. Even the large scales 
are assumed to have moderate gradients only. 
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We show how several stabilized methods enter the framework, both first and 
second order accurate ones. The interest of the first order artificial viscosity method 
is primarily its close relationship to the vertex centered finite volume method. Note 
also that the estimates with an order proposed herein for nonlinear monotone schemes 
to the best of our knowledge are the first of their kind in the literature. 

It appears that for implicit large eddy simulations both the estimate (|1.1[) for 
smooth solutions and the estimate (|1.3[) for rough solutions derived herein are desirable 
properties for the theoretical justification of a method. 

Future work will focus on numerical investigations both in two and three space 
dimensions. Of particular interest is to study the stability of the incompressible Eulcr 
equations to see if the limit estimate with no allowed small scales is sharp. 

Acknowledgment. Partial funding for this research was provided by EPSRC 
(Award number EP/J002313/1). 

Appendix A. Proof of Proposition 12.21 We introduce the discrete errors, 
with Ih denoting the Lagrange interpolant 

eh,i> ■= - and e h ^ '■= u h - it l lu. 

First consider the second equation (|2.4[) and use Galerkin orthogonality 

II Ve fcl *|| a = (V(* - 4*), Ve ft ,*) - (w - u h , e h ^). 

Applying Poincares inequality followed by Cauchy Schwarz inequality we obtain the 
following bound for VP in terms of the error in the vorticity 

l|Ve h ,»|| < ||V(* - h^)\\ + \\u- 7r£W h || + |K.J. 

Consider now the equation (|2.3[) taking Vh = &h,u arL d observing that there holds 

2^ll e 'i,wl| 2 + s i u h; eh,w,eh,w) 

= (d t (uj - 7r L cj), e h , u ) + {u,u - Ve h ^) 

- (tt l uj, u h ■ Ve hjU ) - s(u h ; tt l uj, e h , u )- 

By integration by parts in time we see that the first term on the right hand side is 
zero, by the orthogonality of the L 2 -projection. We then add and subtract Uh in the 
second term on the right hand side to obtain 

+ ttluj, u h ■ We h , u ) - s(u h ; tt l uj, e h)U ) = I + II + HI. 

In the first term on the right hand side we now reintegrate by parts and use Cauchy- 
Schwarz inequality, 

I < \\u\\w l ><*>\\u - «h||||e/,, u || 

< |M| wl ,oo(||V(tf - J h *)|| + \\u - w h ||)||e M || 

< ||o;|| w i,oc(||V(* - I h ^)\\ 2 + n L u\\ 2 + \\e, huJ \\ 2 ). 
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In the second term we use the orthogonality of the L 2 -projcction to retract some 
function Vh and then apply (|2.5[) . 

II = (to - 7TLW, Uh ■ Ve/,., w - Vh) 

< c\\h~i (co - TTLUj)\\s(uh; eh,u>, e/t, w )' 

< ch 2s ~ 1 \\u\\ 2 Hs + -s(u h ;eh,ui,eh,u)- 

For the stabilization term finally we apply the Cauchy-Schwarz inequality and an 
arithmetic-geometric inequality to obtain 

77/ = s(u h ; tt l uj, e h ,u) < s(u h ; it l uj, tt l lo) + -s(uh] e h , u , e h ^). (A.l) 

Then we observe that by adding and subtracting I^rot 'J we may write 

s(u h ; n L u, 7T L w) < s(uh - Jfcrot ttlcj, i: l uj) + s(7fcrot 7r L cj, it l uj) 

and using the definition (|5.1|) and the stability of the L 2 -projection on quasi uniform 
meshes, we have, 

s(u h - I h iot tt l u, tt l uj) < / h 2 \u h - h^ot * | 2 1 V7r L Wft | 2 ds 

< ||Vw|||«, (n) /i||u h - Jfcrot n 2 
< ||V W ||io. C n)MI|V* - 4W)|| 2 + ||V(* - I^)|| 2 + || W - ^ + || e ^|| 2 ) 

and then 

s(I h TOt 7T L W, TT L L0) 

< H4rot nl°°{Q) Y, Hk (l' V ^ _ *^)Wk + h 2 K \\V{u - ^)|| 2 ) 

K 

< \\ u \\'L°°(Q)C'h'K 1 |l a; lli 2 (/:// s (0))- 

We conclude by collecting the upper bounds for the terms I — III, applying approx- 
imability and Gronwall's lemma that 




Here we assumed /i||u;||wi,oo(q) < 1 and neglected the dependence of M| 2 oo(q) (that 
is upper bounded by \\c*}(-, 0)||^„o/q\). It follows that for I = 2 and sufficiently smooth 
solutions we have 

\\(uJ-UJ h )(;T)\\ + \\(u-U h )(.,T)\\<h-z. 
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